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I. INTRODUCTION 

Nuclei near the neutron drip-line provide us with many 
new physics issues which arise from the presence of 
weakly bound neutrons and the coupling to unbound neu- 
tron states. The ground state and the excitation modes 
of a near-drip-line nucleus are indeed very different from 
those of stable nuclei as is testified by the observations 
of the neutron halo[l|, the neutron skinQ and the soft 
dipole excitation Q. In addition the nucleon correlations 
such as the pairing may also be influenced in the new 
circumstances [4j,|5(. Consequently there has been consid- 
erable efforts in the last two decades to develop nuclear 
many-body theories toward this direction. 

Focusing on near-drip-line nuclei in the medium 
mass region, theoretical approaches based on the self- 
consistent mean-field methods or the density functional 
theories are of great promise. The Hartree-Fock- 
Bogoliubov (HFB) theory d , especially those employing 
the coordinate-space representation [J, H 0|, has been 
playing a central role to describe the ground state and 
the pair correlation. The HFB theory provides us also 
with the basis for further theoretical developments to de- 
scribe the dynamics, e.g. the excitation modes built on 
the ground state. Indeed new schemes of the quasiparti- 
cle random phase approximations (QRPA) formulated on 
the basis of the coordinate-space HFB have been recently 
proposed and applied extensively to studies of multipole 
res pon ses of unstable nuclei lp. 13. Hol . [Til . H3 . FPU El EE 
EE US El El, HE HE HI (see also references 

in Ref.H). 

There are two important requirements to be considered 
when the HFB+QRPA theories are applied to near-drip- 
line nuclei. First of all, the coupling of excitation modes 
to the continuum states have to be taken into account 
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since most of excitation modes including even the low- 
lying excitations are located near or above the nucleon 
separation energy. This can be achieved by means of the 
continuum QRPA methods 0, EH, El El El, El • Sec- 
ondly, the QRPA description should be consistent with 
the HFB description of the ground state in the sense 
that the same effective interaction or the same density 
functional should be used for both descriptions. If this 
is achieved, one can calculate the ground and excited 
states solely from the effective interaction (or the energy 
density functional) without relying on phenomenologi- 
cal parameterization of the mean-fields. This is often 
called the requirement of the self-consistency. The two 
requirements, however, have been in a trade-off relation 
in the actual implementations. Namely in the contin- 
uum QRPA methods which fulfill the first requirement 
the self-consistency has been left behind since the resid- 
ual interaction for the QRPA description is often approx- 
imated to a tractable simple contact force [H flfl. Il2l| or 
the Landau- Migdal forces [H El El- On the other 
hand, recently developed fully self-consistent QRPA's 
using the Skyrme function al up . l24j and the relativis- 
tic mean- field functional [let, Il9l| treat approximately the 
continuum states by employing the finite-box discretiza- 
tion or the discrete oscillator basis. 

It is therefore important to develop a new formulation 
of the continuum QRPA which is based on the nuclear 
density functional and thus satisfies the self-consistency 
as precisely as possible. In the present paper we try to 
make a one step progress in this direction. 

To this end we shall proceed in the following way. 
We start with the Skyrme's Hartree-Fock energy func- 
tional combined with the pair correlation energy. We 
then use this functional not only for the static Hartree- 
Fock-Bogoliubov mean-fields but also to derive the resid- 
ual interaction to be used in the continuum QRPA. In 
formulating the new continuum QRPA, we pay special 
attention to the energy weighted sum rule, which is not 
satisfied in the previous continuum QRPA's [IE EH, El, 
El El El ■ To satisfy this we take into account explicitly 
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the velocity dependent central terms of the Skyrme ef- 
fective interaction to derive the residual interaction, and 
then implement the residual interaction into the Green's 
function formulation of the continuum QRPA proposed 
in Ref. 10]. In the present paper, however, we do not in- 
clude the spin dependent densities and the Coulomb in- 
teraction in deriving the residual interaction, and hence 
the goal of the full self-consistency is not achieved yet. 
Our formulation of the Skyrme QRPA is similar to that 
of Ref. [20L |2~H apart from the treatments of the contin- 
uum quasiparticle states, on which we impose the out- 
going wave boundary condition instead of the finite-box 
discretization. 

By performing numerical calculations, we shall demon- 
strate that the Skyrme continuum QRPA in the present 
formulation indeed satisfies the sum rule as far as the 
dipole and quadrupole responses with natural parities 
are concerned. We shall also show that the inclusion 
of the velocity dependent terms gives better description 
of the strength function and the transition densities of 
the multipole responses in comparison with the previ- 
ous continuum QRPA that utilizes residual interactions 
of the simple contact forces. 

The paper is organized as follows. The formulation is 
given in the next section. In Section III we present re- 
sults of numerical calculation performed for the isovector 
dipole and the isoscalar/isovector quadrupole responses 
in neutron-rich O and Ca isotopes. We discuss both 
the low-lying excitations and the giant resonances. We 
shall illustrate in detail the importance of taking account 
of the velocity dependent terms by comparing with the 
Landau-Migdal approximation of the Skyrme effective 
interaction[27l [28|. The conclusions are drawn in Sec- 
tion IV. 



II. CONTINUUM QRPA USING THE SKYRME 
FUNCTIONAL 



In this section, we give a formulation of the continuum 
QRPA which is based on the Skyrme functional. 

We start with the energy functional of the system de- 
fined for a determinantal many-body state vector 
of the generalized form in which the pair correlation is 
taken into account by means of the Bogoliubov's quasi- 
particle method The time-dependence is explicitly 
written here since we consider dynamical multipole re- 
sponses of the system under a time-dependent perturba- 
tion. The energy functional E — Eskyrme+E pa i r consists 
of the Skyrme Hartree-Fock energy Eskyrme and the pair 
correlation energy E pa i r . The Skyrme Hartree-Fock en- 
ergy Eskyrme is expressed in terms of the local density 
p q (r,t), its spatial derivatives V/9 g (r,i) and Ap q (r,t), 
the current density j q (r,t), the kinetic energy density 
T q (r,t), the spin density s q (r,t) and the spin-orbit ten- 
sor J q (r,t) where q = n,p stands for the neutron or 
proton components p9l. [30j . Given a parameter set such 
as SIIl[3ljl, SkM*[lJand SLy4[H, the Skyrme Hartree- 



Fock energy functional E S k y rme[p, Vp, Ap, r, j, s, J] is 
completely specified. Concerning the pair correlation en- 
ergy Ep a i r , we use the one evaluated for the the density- 
dependent delta interaction (DDDI) (34|, l35j 



V pair (l,2) = -V (l-P„) 



>1 



p(r) 



5{r-r').{\) 



Epair is a functional of the local density p q (r,t) and the 
local pair densities 

P± q (r, t) = (*(t)| ^(r |)V4(r ])±^ q (r ^ q (r j) |$(t)) . 

Application of the static variational principle to the 
total energy functional Eskyrme 
Hartree-Fock-Bogoliubov equation 



Ep a i r leads to the 



W O9 ? (rcr) = E q <p q (ra) 
for the quasiparticle wave function 



(j> q {ro) 



Here 



7~taq — 



<Pq,l ( rCr ) 
<Pq,2 {TO) 



-h* q + X q 



(3) 



(4) 



(5) 



is the 2x2 matrix representation of the HFB mean-field 
Hamiltonian 

h = I drdr'^2h q (ra,r'a')i()l(ra)i() q (r'<T') 



1 

+ 2 
+h.c 



drdr' h q (ra, r' a')ip^ q (ra)ipl('r'' <r') 



(6) 



The Hartree-Fock Hamiltonian h q and the pair potential 
h q are defined through the functional derivative of the 
energy functionals Eskyrme and E pa i r , respectively. 

We consider multipole response of the nucleus under a 
small time-dependent external perturbation 

V ext (t) = e-**Y, j drf q (r)J2^q(.ra)4, q (r<j) 



-h.c. 



(J) 



expressed in terms of a one-body spin-independent local 
field fq(r), for which we take a multipole field cx t l Ylm 
such as the electric dipole and the isoscalar/isovector 
quadrupole fields. 

The external perturbation causes the induced fields in 
the Hartree-Fock mean-field and the pair potential, which 
we denote 5h q and Sh q , respectively. Sh q and 5h q are 
expressed in terms of fluctuations in the various one-body 
densities 

Sp q (r, t), SVp q (r, t), SAp q (r, t), Sr q (r, t), Sj q (r, t), 
SJ q (r,t),Ss q (r,t), 

Sp± q (r,t) (8) 



3 



and the second derivatives of the energy functional. A 
fully self-consistent QRPA based on the Skyrme HFB 
functional can be constructed if one considers all the 
kinds of density fluctuations in Eq.®. In the previous 
continuum QRPA approaches, however, only the fluctu- 
ations in the local densities Sp q (r,t) and 5p± q (r,t), and 
the induced fields associated with these density fluctua- 
tions have been taken into account [ToL [TTl [T^ [T^. [l5l| . 
Although this approximation has a large number of prac- 
tical usefulness, it is not sufficient in some respects: it 
violates the energy weighted sum rule when the Skyrme 
HF mean-field with the effective mass is adopted. This 
is because the current conservation law is not satisfied 
when the velocity dependent parts (the terms propor- 
tional to the t\ and ti terms) of the Skyrme interac- 
tion and the current fluctuations Sj q are neglected in the 
RPAm^Hii]. We aim at improving this point. For 
this purpose we shall include all the density fluctuations 
that are responsible for the energy weighted sum rule for 
the responses caused by the spin independent local mul- 
tipole fields. They are the fluctuations 5j , SVp q , SAp q , 
and Sr q in the current, the spatial derivatives of the den- 
sity and the kinetic energy density. It is more preferable 
to take into account also the fluctuations in the spin- 
dependent densities s q and J q , but we neglect them in 
the present work. This is one approximation which re- 
mains in the present approach. We neglect the residual 
Coulomb interaction. The cross derivatives among p q and 
p± q are also neglected. These are the second approxima- 
tion we introduce. Consequently the full self-consistency 
is not fulfilled, but the treatment of the residual interac- 
tion is significantly improved com par ed to the previous 
continuum QRPA approaches [H [ill, El El Ei El in 
the sense that the present formalism allows us to describe 
the correct energy weighted sum rule for the multipole 
responses. Note also that the approximate treatment 
of the particle-hole residual interaction is comparable to 
that adopted in the currently available continuum RPA 
approaches which utilize the Skyrme-Hartree-Fock func- 
tional without taking into the pairing [H EH El EH- 
On the other hand, we should keep in mind that the ap- 
proximation neglecting the spin dependent densities may 
not be justified for multipole responses with unnatural 
parities involving spin excitations. 

The induced fields under the above approximations are 
expressed as 



Sh q = a qq i5p q > + b qq >5A + p q , 



A + A 



b qq '5p q ' 



,8V Pq , 



■b qq ,62ij q , (9) 



and 



Sh q = a q 5p +q - a q Sp 



(10) 



The functions a qq > , b qq > , c qq > and a q are expressed in terms 
of the local densities and the effective interaction pa- 
rameters. The definition of a qq >,b qq > and c qq > follows 
Ref. [43l | , and their detailed expressions are given in Ap- 
pendix. Note here that we use SA + p q (r,t) = (A + 
A')8p q {rr',t)\ r , =r = 5Ap q (r,t) - 25r q (r,t) in place of 
5Ap q (r,t), where p q (rr',t) is the density matrix. We 
also attached a factor 2i to the current fluctuation Sj q so 
that S2ij q = (V — V')8p q (rr' , t)\ 1 becomes in par- 
allel with 5Vp q = (V + V')8p q {rr'~t)\ r , =r . The fluc- 
tuation in the kinetic energy density Sr q (r,t) does not 
appear here since it can be eliminated from the induced 
field by a partial integration. Nevertheless we eventually 
include Sr q as explained later. 

The induced fields are also represented in the 2x2 
matrix form as 



Sh q Sh q 
5K -Shi 



(11) 



where Sp y is a collective notation for 



Sp 1 G 8pq.5A + p q ,8V p qi 82ij ql 5p±q. 



(12) 



and O 13 denotes the derivative operators 1, A~+ A*, ^ — ^ 
and ^ + ^ while B@ stands for one of the 2x2 matrices 

1 Wi 0W0 i\ , ( 1 



-1 



1 



1 



and 



-1 



np a rep- 



resents the functions a qq ' , b qq i , c qq 



and a q in Eqs.([9]) and 
(fTU)) . The correspondence among O a , B a np^ and Sp 7 is 
shown in Table [J 

The external perturbation and the induced fields cause 
quasi-particle excitations, which in turn bring about 
fluctuations in the densities Sp qi Sp± q ,8A + p q ,V8p q and 
S2ij q . This relation is given by the linear response equa- 
tion, which is written as 



5p a (r,u>) = J2 J dr'R^(rr\ 



(13) 



in the frequency domain. 



Here RQ q (rr',uj) is the unperturbed response function for the density Sp a and the field B^O 13 . Using the Green's 



4 



function formalism of the continuum QRPA[10(, the un- perturbed response function is expressed as 
I 



47TI 



+ 



47TZ 



dETr \A a d a {r)g 0q (rr', E + %w + ie)B O p (r')g 0q {r'r, E) 
dETr \A a d a (r)g 0q (rr', E)B f) 6 f3 (r l )g 0q {r'r : E-Hu- ie) 



(14) 



in terms of the quasi-particle Green's function g^ q {E) = 
(E — Hoq)^ 1 and a contour integral in the complex en- 
ergy plane. The complex energy integral is performed on 
a rectangular contour C enclosing the negative energy 
par t of the real E axis with the two sides located at ±|e 
hfl. Here e is a small parameter which plays a role of 
the smoothing energy width. The matrices A a and B 13 
and the operators O a and O 13 follow Table U but we re- 
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mark that the matrix A a takes a form A a 







for 



the particle-hole densities 5p q ,8A+p q ,SWp q ,52ij q and 
8r q while the matrix B" has the following definitions: 

q ^ ~\ for the 'time-even' quantities Sp q , 5A + p q 

and S'Vpq, and B 13 = f q ^ ^ for the 'time-odd' <52zj g 
(See Table I|. 

Let us assume the spherical symmetry of the ground 
state, and we apply the multipole decompositions. Let 
L be the multipolarity of the excitation modes under 
consideration. The fluctuation in the scalar quantities 
S p a = Sp q , Sp± q and SA + p q are expanded as 



Sp a (r,u;) = Y LM (r)[Sp a } L /r 2 



(15) 



and we now consider only the radial functions [J/j q ]l — 
[5p q ]L, [6p± q ]z, and [5A+p q ]L. Concerning the vector 



quantities 5p a = <5V/9 9 and 52ij q , they are expanded 
as 

6p a (r,u)= Y LXM (r)[5p a } x L /r 2 } (16) 

A=L±1 

in terms of the vector spherical harmonics Y lxm 
and the radial functions [<5p Q ]£ = [S'V p q ]^~ L±1 and 
[S2ij q ]l =L±1 . Note that only the terms A = L ± 1 re- 
main here since we consider the multipole excitations 
with the natural parity. Then the linear response equa- 
tion (|13p is rewritten as an equation for the relevant 
radial functions [Sp q ] L ,[Sp± q } Ll [SA + p q } L ,[5Vp q }l =L±1 
and [82ij x . Denoting collectively these density 

fluctuations Sp a L, the linear response equation for these 
variables is given by 



5p aL (r,u>) 

= J2jdr'Rf qL (rr'c) 



^2K 01 8p lL (r',uj)/r /2 + 5p fiq f qL {r') 



(17) 

using the unperturbed response function for the fixed 
multipolarity L 



Airi l~ dE ^ 



-Tr 



Ijl'j' 



\(i'j'\\Yl\\Ij)\ : 

2L + 1 



Tr 



A a 6f oVr {r)g . qVf {rr',E + hu; + i^O^^Q^ir'r, E) 



A a Op jnj (r)g 0M (rr',E)B^6f jl , j , {r>)g , qVj , (r'r, E - hu; - ie) 



(18) 



Here go,qij{r'r, E) is the 2x2 radial HFB Green's func- 
tion for specified orbital and total angular momenta I 
and j, and 0^ .,y is the radial derivative operator cor- 



responding to the previously defined O". Their explicit 
forms are given in Table HI We adopt the exact form for 
the radial HFB Green's function! 441 constructed as 



5 



G , qlj (rr',E) = £ c^{E) (o(r 



(19) 



r 



in terms of two independent solutions <p^ s ) (r, E)(s = 1,2) 
regular at the origin r = of the radial HFB equation 
and two independent solutions (^y (r, E)(s — 1,2) satis- 
fying the out-going boundary condition. (Th e construc- 
tion (fT9|) is the same as that used in Refs.[l(l 13 except 
that the effective mass should be taken into account in 
the definitions of the Wronskian and the coefficients 
while in Rcfs. 

0,13 

the bare mass is assumed.). In this 
way the exact treatment of the continuum single-particle 
states satisfying the proper boundary condition is imple- 
mented in the QRPA formalism. 

Note that in Table Q] we use the following convention 
for the derivative operators marked with the right/left- 

When this derivative 
of 

~8 



n n 

sided arrows such as -fe- ± -fe- 
ar or 

is inserted in 0^ jV y(r') in the first term of r.h.s 



Eq. (|T8]) . the derivative symbol -*jp with the right-sided 
arrow indicates that it acts on the coordinate r' in 
the Green's function Go,qij (r'r, E) while the other one 

4rj with the left-sided arrow acts on the Green's func- 

or' 

tion Go,qi'j'(rr',E + Tiuj + ie). The same rule is ap- 

plied also to the operator i-e., j^r acts on r 

7T 

in Go, q i> j {rr' , E + hu; + ie) while on r in Go, q ij (rr, E) . 
The ordering of the operators and the Green's functions 



makes sense in Eq. (fTS| . 

To obtain a numerical solution of the linear response 
equation, we need to rewrite further Eq.(fl7|). When the 
radial derivative operators and -gp- act on the radial 

HFB Green's function like -^pGo, q ij(rr' , E)-gp, a singular 



2m* ( r) 

term proportional to — ^ — S(r 



emerges. We need 



to treat these singular terms separately in the numerical 
calculation. For this purpose we rewrite the derivative of 
the Green's function into singular and regular parts 



2m* q (r) 



—p ! -S{r-r r ) 
n 

~Go, q ij(rr',E)-^ (20) 



where the regular part (the second term in r.h.s de- 
noted with the tildered derivatives) is defined as a part 

of -§pGo,qij (rr' , E)Jpp that arises from the action of 



{+*) 
ij *" ,M vij 
but not on the Heaviside theta function 



and 



the derivatives on the wave functions 
in Eq.([T£ 

9(r — r'). Inserting this decomposition into the response 
function(Eq. (fl8|) ) . the r.h.s. of the linear response equa- 
tion (fTT|) is decomposed into two parts: 



J 



5p aL {ruj) =^ Jdr'R* P qL {rr'uj) 



«/3 7 {r')Sp jL (r'uj)/r' 2 + Sp. 0q f qL (/) 



+2^Sf (r) 





^K / 3 7 (r)<5p 7L (rcj)/r 2 + S 0tOq f qL (r) 



(21) 



where 



S PaL G [Spqh, [SA +Pq ] L , [5V Pq ] x = L±1 , [Sr q ] L , [82ij q f= L±l , [6p ±q ] L . 



(22) 



Here Rq^l denotes a part of the response function 
which contains only the regular parts of the derivatives of 
Go.qij ■ Its expression is the same as that of Rq qL (Eq. (fT8]) ) 
except that the derivatives of the Green's function, e.g. 

■^Go, q i 3 (rr', E)-Qj in Eq.fgOl) is replaced by the corre- 



sponding regular part j^Go^ij (rr' , E)-Qj 



hand, the second term of r.h.s. of Eq. pT 



. On the other- 
represents con- 



tribution from the singular terms such as — p — S(r- 



<-')■ 



The integral J dr' disappears in this term because of the 
delta function. The expressions of S q @ are given in Ap- 
pendix. Note that S q @ is a one-point function indepen- 
dent of the frequency uj, expressed in terms of local quan- 
tities such as p q {r), T q (r), m*(r) and their derivatives. 

It is noted that the linear response equation (|2"Tj) in- 
cludes the fluctuation [5r q ]£ in the kinetic energy density 



6 



6 p a 



O a {r) 



5p q 



5A+p q 



SVp q 



2 



' -l 



1 



A + A 



dr' 2 dr' 2 



-pr 



l i d , a i L-i 

2L+1 Sr + 9r ' r 



(for [*Vp,]*= 



<5t„ 



2 




1 
1 



d d i _d Q_ i 

dr dr r dr dr r 

| Ui+l)+i'(i' + l)-(L+2)(L-l) [J_ 



V^(I-I+ ' ( ' +1) '/ ( '' +1) f) (far laayj^" 1 ) 



_L 

2L+ 



<5p 4 



1 

1 



Q>q3qq' 



Sp^ 



1 
-1 



TABLE I: The correspondence and the expressions for the matrices A a and S , the operators 0° 
in Eqs. (fTTj) . (fl4|) and (fl8)l . See also the text and Appendix. 



and d^j'j' and K a/ s appearing 



T q as a dynamical variable to be considered. This is be- 
cause [<$Tq].L emerges from the singular terms associated 
with the linear response equation for [5A+Pq]i,. Finally 
we make a little remark on the structure of the singu- 
lar terms. The presence of the singular terms has been 
notified in the formulation of the Skyrme-HF plus contin- 
uum RPAptol, |43| where the pairing is neglected. In the 
present Skyrme-HFB plus continuum QRPA approach, 
the structure of the singular terms is more involved since 
the response function contains two single-particle HFB 



Green's functions (instead of one Green's function in the 
case of the continuum RPA[4(| EH). Looking at the ex- 
pression of Eq. (fT5)l , it may appear that products of two 



delta functions 



2m; (r) 



5(r — r') emerge from the singular 



terms of two HFB Green's functions in Eq. (fl8|) . Such a 
term however does not contribute to the response func- 
tion since it has no energy dependence and hence it van- 
ishes when the contour integral in the complex energy 
plane is performed. 



III. NUMERICAL ANALYSIS 



A. Numerical procedure 



In this section, we shall demonstrate the Skyrme con- 
tinuum QRPA by performing numerical calculations for 
the dipole and quadrupole responses in 20 O and 54 Ca. 



Let us first describe the detailed procedure of the nu- 
merical calculation. 

We adopt the SkM* parameter set of the Skyrme 
interaction and the mixed-type parametrization of the 
DDDI pairing interaction (rj = 0.5, 7 = 1, po = 0.16 
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fm -3 ) 35] for most of the calculations. The force strength 
Vq of the DDDI is chosen so that the average neutron 
pairing gap (A n ) reproduces the overall magnitudes of 
the experimental odd-even mass differences for the iso- 
topic chain, obtained with the three-point formula (45|. 
Here we use the average pairing gap defined by (A n ) = 
J drp n (r)A n (r) / J drp„(r). The adopted value is Vq — 
-280 and -285 MeVfrn" 3 for 20 O and 54 Ca producing 
(A n ) = 1.91 MeV and 1.29 MeV, respectively. 

Since we use the contact interaction for the effective 
pairing interaction, we need a cut-off of the quasi-particle 
states in the HFB calculation. We define the cut-off with 
respect to the quasi-particle energy E a < E max = 60 
MeV. Concerning the angular momentum quantum 
numbers Ij we sum up the quasi-particle states up 
to Imax — 7h and 8h for 20 O and 54 Ca, respectively. 
In performing the HFB and the continuum QRPA 
calculations, we discretize the radial coordinate space up 
to r max = 15 fm with an equidistant interval Ar = 0.2 
fm. In the continuum QRPA calculations, the dynamical 
quantities to be obtained are the eighteen functions 
[Sp q ] L , [Sp± q ] L , [SA +Pq ] L: [8V Pq \ x = L± \ [52ij q ] x = L±l 
and [ST q ]L which obey the linear response equation 
([2"Tj) . Using the same radial mesh, these functions 
are represented as a grand vector while the linear 
response equation is represented as a linear algebraic 
equation where the response function Rq qL (and S q @) 
corresponds to a matrix. Since the number of the 
functions to be solved is larger (18 vs. 6) than in the 
previous continuum QRPA that handles only the local 
densities and [Sp± q ]L, the number of the matrix 

elements of the response functions is therefore about 
ten times larger than in the previous continuum QRPA 
calculations. To reduce the increased computational cost 
thus caused, we have chosen the values of l m ax and r max 
smaller than those used in our previous calculations 
[ToL [Til . fl2L [HI . For the same reason we have used here 
a relatively large smearing parameter e = 1.0 MeV in 
most of the following calculations. We evaluate the 
strength function at discretized excitation energies with 
an interval of 0.5MeV. 

It is noted here that the self-consistency is not com- 
pletely satisfied in the present formulation since a few 
approximations are introduced in deriving the residual 
interaction from the Skyrme HFB functional. Conse- 
quently the spurious modes of motion which should have 
exact zero excitation energy according to the Thouless's 
theorem^ do not emerge at the expected energy. A com- 
monly adopted procedure to circumvent this problem is 
to renormalize the residual interaction n a p in Eq. (|21|) 
by an overall factor / as K a p — > f x n a p so that the 
excitation ener gy o f the spurious mode is forced at the 
zero energy[IE III El El EI M, M, M, El . We ap- 
ply this renormalization procedure to the particle-hole 
residual interactions that are derived from the Skyrme 
HF functional Eskyrme- The residual interaction in the 
particle-particle channel, derived from the pair correla- 
tion energy E pa i ri is kept in the original strength since 
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FIG. 1: The B(E1) strength function of isovector dipole re- 
sponse in 20 O (upper panel) and in 54 Ca (lower panel) cal- 
culated with the parameter set SkM*. The solid curve is the 
result obtained in the full calculation while the dashed curve is 
that in the Landau-Migdal (LM) approximation. The dotted 
curve in the upper panel is the one in the to + tz approxima- 
tion. See also the text. 



the continuum QRPA in the Green's function formalism 
fulfills the self-consistency in the particle-particle channel 
with high accuracy [lol EI| ■ The renormalization factor 
is / = 1.0470 and 1.0142 for 20 O and 54 Ca, respectively. 

In the following analysis, we would like to demon- 
strate how the description of the multipolc response 
is improved in comparison with the previous contin- 
uum QRPA where the residual interaction is simpli- 
fied to a contact force. For this purpose, we per- 
form calculations where the Landau-Migdal (LM) ap- 
proximation to the residual interaction is introduced 
IH El, EE EE E3, HI U HI- This is an approxima- 
tion which replaces the residual interaction by a contact 
force oc S(r — r') whose strength is given by the density- 
dependent Landau-Migdal parameters Fq and Fq evalu- 
ated for the Skyrme functional [E HE [28[ using the local 
density approximation. It should be noted, however, that 
the Landau-Migdal parameters Fq and Fq contain a part 
of the t\ and ti terms, and this approximation should be 
distinguished from dropping all the t\ and £2 terms. 
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FIG. 2: The same as FigQ] but for the B(IS2) isoscalar 
quadrupole strength function in 20 O and 54 Ca. 
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FIG. 3: The same as Fig[TJ but for the B(IV2) isovector 
quadrupole strength function. 



B. Strength function 

The strength function 
S(hu) = Y,\( v \^M\0)\ 2 S(huj-E,) 

v,M 



—^— Im Yl J dr fiL(r)[Sp q ] L 



(r )W )(23) 



for the operator Flm with the multipolarity L can be 
evaluated in terms of the solution [5p q ]L(r,u>) of the 
linear response equation (|21[) obtained for the external 
field F L0 - We evaluate the S(E1), B(IS2) and B{IV2) 
strength functions associated with the electric dipole op- 
erator 



ftiv 



eN Z eZ N 

—j- nYlM(Cli) -r J]] TiYiM^i 

i=l i=l 



(24) 



and the isoscalar/isovector quadrupole operators 



pis 

r 2M 



(25) 



The B(E1) strength functions calculated for 20 O and 
54 Ca are shown with the solid curve in FigfTJ The broad 
peaks around E x = 20 MeV in 20 O and E x = 16 MeV 
in 54 Ca correspond to the giant dipole resonance (GDR). 
There is s small bump around E x — 8 MeV in 54 Ca, 
which corresponds to the soft dipole excitation or the 
pygmy dipole resonance. We find however that the small 
peak E x = 13 MeV in 20 O is neither the GDR nor the soft 
dipole excitation, but rather a non-collective two quasi- 
particle excitation (cf. Section Till D[) . The strength at 
E « is due to the spurious mode, and it is caused by 
the incomplete self-consistency. 

For the sake of comparison, the B(El) strength func- 
tion obtained in the Landau-Migdal (LM) approximation 
of the residual interaction is also plotted with the dashed 
curve in FigfT] Note that the renormalization factor used 
in the Landau-Migdal approximation (/ = 0.6686 and 
0.7515 for 20 O and 54 Ca, respectively) deviate signifi- 
cantly from one in contrast to those in the full calculation 
(/ = 1.0470 and 1.0142). This fact suggests that there is 
significant improvement in the self-consistency compared 
with the LM approximation. This feature is pointed out 
in a Skyrme-QRPA calculation using discretized contin- 
uum quasiparticle states [21]. 

It is seen in FigfT] that the profile of the strength func- 
tion obtained in the LM approximation differs signifi- 
cantly from that in the full calculation. The peak po- 
sitions of the giant dipole resonance are apparently dif- 
ferent. Estimating the centroid energy of the GDR by 
-B(GDR) = mi/mo using the energy weighted sum mi 
and the non-weighted sum mo, we find F(GDR) = 20.66 
MeV ( 20 O) and 15.80 MeV ( 54 Ca) for the full calculation 
while F(GDR) = 17.67 and 14.20 MeV in the LM ap- 
proximation, exhibiting a rather large difference by about 
2-3 MeV. It is clear that the LM approximation is not 
very appropriate to give precise quantitative description 
of the GDR. 

If we evaluate the energy weighted sum integrated up 
to E = 15 MeV for 20 O, the full calculation gives 9.7 % of 
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the classical Thomas- Reiche-Kuhn (TRK) sum rule value 
while it is 20.1% in the LM approximation. Comparing 
with the experimental value 12% we find that the full 
calculation is in better agreement with the experiment. 
The B(E1) strength function in 20 O is calculated in a 
fully self-consistent Skyrme-QRPA calculation 24] using 
the same SkM*. We find only small difference between 
our calculation and that in Ref . [24| . It may be attributed 
to the neglect of the Coulomb and spin-dependent terms 
in our calculations. The observed effect of the velocity 
dependent terms on the GDR centroid energy is essen- 
tially the same as that discussed in Ref. [2lT] . 

Figure [2] displays the £?(IS2) isoscalar strength func- 
tion for the quadrupole responses in 20 O and 54 Ca. In 
both nuclei there are two significant peaks, one around 
E x = 2 — 3 MeV corresponding to the low- lying 2 + collec- 
tive vibrational mode and the other around E x = 15 — 20 
MeV corresponding to the isoscalar giant quadrupole res- 
onance (ISGQR). (The experimental 2f energy in 20 O 
is 1670 keV|48j.) The calculated isoscalar quadrupole 
strength function for 54 Ca is quite similar to that ob- 
tained in the fully self-consistent Skyrme QRPA using 
the same SkM*[2J] apart from features associated with 
different choices of the smoothing width. 

Concerning the effect of the velocity dependent terms, 
it is seen in Fig.[2]that the difference between the full cal- 
culation and the LM approximation is less significant in 
comparison with the isovector dipole response: the peak 
positions of the giant isoscalar quadrupole resonance 
(E = 19.0 MeV) and of the low-lying state (E = 3.0 
MeV) in 20 O is affected only little by inclusion of the ve- 
locity dependent terms. The same is seen also in 54 Ca. 
Note however that the influence of the velocity depen- 
dent terms on the £?(IV2) isovector distribution is clearly 
larger than in the case of the -B(IS2) isoscalar strength 
distribution as is seen in Fig[3] Combining Figs. [TJ [2] 
and [31 we see an apparent trend that the influence of 
the velocity-dependent terms is more significant in the 
isovector responses than in the isoscalar responses. 

Before moving to the next subsection, we would like to 
make a few additional remarks on the effect of the veloc- 
ity dependent terms. We first remark that it is possible to 
consider another way to evaluate the effect of the velocity 
dependent terms, e.g. by comparing with a calculation 
where all the velocity dependent terms containing the t\ 
and t2 parameters are completely neglected. (In other 
words, it is the calculation where only the simple contact 
interaction associated with the to and £3 terms are taken 
into account. It is different from the Landau-Migdal 
(LM) approximation since in the latter a part of the t\ 
and t% terms is renormalized into the Landau-Migdal pa- 
rameters Fq and Fq.) The -B(El) strength function in 
this to + ts approximation calculated for 20 O is plotted 
in the upper panel of Fig[T] together with the other two 
curves representing the full calculation and the LM ap- 
proximation. We find here that the result obtained in the 
to + £3 approximation is almost identical to that in the 
LM approximation. Secondly, we remark that the effect 
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FIG. 4: The same as the upper panel of Fig[TJ but the Skyrme 
parameter set SLy4 is used. 



of the velocity dependent terms depends on the adopted 
Skyrme parameter set. To demonstrate this we show in 
Fig. 4 the B(E1) strength function in 20 O obtained with 
SLy4[l| instead of SkM*. It is seen that there is no big 
difference in the GDR peak position between the full cal- 
culation and in the Landau-Migdal approximation, and 
hence the effect of the velocity dependent terms in the 
case of SLy4 appears smaller than in the case of SkM*. 
Note however that even in this case there is significant 
difference between the LM and to +t^ approximations. If 
we look at the difference between the full calculation and 
the to +^3 approximation, the effect of the velocity depen- 
dent terms is not negligible. Note also that the difference 
between the full calculation and the LM approximation 
is not negligible in the sum rule (cf. next subsection). 



C. Energy weighted sum rule 

Let us analyze whether the energy weighted sum rule 
is satisfied in the present calculations. For this purpose 
we evaluate the running energy weighted sum defined by 



W{E X ) 



dE E S(E), 



(26) 



which integrates the sum up to an excitation energy E x . 
The limiting value of W(E X ) for a sufficiently large E x 
is to be compared with the energy-weighted sum rule 
(EWSR). 

The EWSR for the -B(ISL) isoscalar multipole strength 
function is identical to the classical sum rule 



,d — 



4tt 2m v 



l) 2 A{r 



2L-2\ 



which is expressed in terms of 



the expectation value of the radial moment r 



2L-2 



with 



respect to the ground state [g, [38|. For the isovector mul- 
tipole strength functions, however, the EWSR contains 
the enhancement factor which arises from the residual in- 
teraction, in particular the velocity-dependent terms in 
the case of the Skyrme effective forcc[36, 38, 39, 49, 50]. 
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The EWSR for the £?(E1) electric dipole strength func- 
tion is given by (El) = ^-mf (1 + k) where K is 
the enhancement factor which is easily evaluated in the 
case of the Skyrme force [H, El- The value of K for SkM* 
is k = 0.32 and 0.36 in 20 O and 54 Ca, respectively. 

The upper panel of the Figj5] displays the running en- 
ergy weighted sum W(E X ) for the B(E1) strength func- 
tion of the dipole response in 20 O (cf. FigQ]). The 
running sum evaluated at the highest calculated energy 
E x = 55 MeV reaches 96% of the EWSR. This suggests 
that the EWSR is satisfied in the present calculation. It 
is noted that W(E X ) approaches the EWSR value includ- 
ing the enhancement factor, but not the classical TRK 
value plotted with the dotted horizontal line in Fig. [5] 
Namely the effect of the velocity dependent terms in the 
residual interaction is indeed included in the present cal- 
culation. In the same figure, we also show the result ob- 
tained in the LM approximation. In this case, however, 
the running sum reaches only 86 % of the EWSR, and 
hence the approximation fails to describe the EWSR and 
the enhancement factor. The fluctuations SA+p, JVp, 
St and, in particular, S2ij play the essential role to re- 
store the EWSR since these are the fluctuations associ- 
ated with the velocity-dependent terms (oc t x and t 2 )- 
Note that the LM approximation neglects these fluctu- 
ations although a part of the velocity-dependent t\ + t% 
terms is taken into account via the Landau-Migdal pa- 
rameters Fq and Fq. 

The lower panel in Fig[5] displays the running en- 
ergy weighted sum for the £?(IS2) isoscalar quadrupole 
strength function in 54 Ca (cf. Fig|2]). The energy 
weighted sum amounts to 98% of the EWSR, and we 
confirm more clearly than in the analysis of the elec- 
tric dipole strength that the EWSR is satisfied in the 
present calculation. When we adopt the LM approxi- 
mation where the fluctuations 5A + p, 5Vp, St and S2ij 
are neglected, the running sum W(E X ) at the maximum 
energy E x = 55 MeV overshoots the EWSR by about 
ten percent. Again the consistent inclusion of the veloc- 
ity dependent terms in the QRPA description is essential 
to guarantee the EWSR. It is noted also that the differ- 
ence between the full calculation and the Landau-Migdal 
approximation is significant mainly in the energy region 
(E x > 18 MeV) higher than the giant resonance peak. 

We confirmed that the EWSR is satisfied also in the 
case of the other Skyrme parameter set SLy4. For the 
£?(E1) and _B(IS2) strength functions in 20 O, the running 
sum W(E X ) at the highest energy E x = 55 MeV reaches 
95% and 97% of the EWSR, respectively. In the LM 
approximation, on the contrary, the sum overshoots the 
EWSR by 6% and 12% in the B(E1) and BQS2) strength 
functions, respectively. 

In FigJHl we demonstrate influence of the residual pair- 
ing interaction on the energy-weighted sum. When the 
residual pairing interaction is dropped in the continuum 
QRPA calculation (i.e., the dynamical pairing effect is 
neglected), the energy- weighted sum overestimates the 
EWSR value by 5%. As pointed out already [H EL 01, 
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FIG. 5: The running energy weighted sum for the £>(E1) 
electric dipole strength function in 20 O (upper panel) and 
that for the _B(IS2) isoscalar quadrupole strength function in 
54 Ca(lower panel) obtained using the Skyrme parameter set 
SkM*. The solid curve represents the result of the full calcula- 
tion of the present Skyrme continuum QRPA while the dashed 
curve is the result obtained in the Landau-Migdal (LM) ap- 
proximation to the residual interaction. The horizontal solid 
line indicates the value of the energy-weighted sum rule in the 
present calculation, which includes the enhancement factor in 
the isovector response. The dotted horizontal line in the up- 
per panel denotes the value of the Thomas-Reiche-Kuhn sum 
rule. 



inclusion of the dynamical pairing effect is important to 
guarantee the energy weighted sum-rule because other- 
wise the self-consistency in the pairing channel would be 
violated. This means, in the present context, that we 
need to include both the velocity dependent terms of the 
Skyrme effective interaction and the residual pairing in- 
teraction. 

Figures [7] and [5] show the B(E1) and 5(IS2) strength 
functions in 20 O, respectively, calculated with use of a 
small smearing width e = 0.2 MeV. They are compared 
with those obtained with e = 1.0 MeV (cf. Figs [1] and 
[2]). The running energy weighted sum of the strength 
function is also shown. Finer structures of the strength 
functions are visible here: we can distinguish the first 
excited 2 + , which is a bound discrete state located below 
the separation energy. 

It is seen from the bottom panels of Figs. [7] and [8] 
that the agreement with the energy weighted sum rule is 
improved with use of the smaller smearing width. The 
running sums at the highest energy E. x = 55(50) MeV 
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dipole response of 20 O using SkM*. The solid and dashed 
curves are the results obtained in the full calculation and 
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dotted curve is the one in which we neglect the residual pair- 
ing interaction in the continuum QRPA. 
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are 98 and 98 % for the B(E1) and 5(IS2) strengths, 
respectively, which should be compared with the corre- 
sponding values 96 % and 96 % obtained with e = 1.0 
MeV. 

In the following we return to e = 1.0 MeV since calcu- 
lations with e = 0.2 MeV demand very long computation 
time. 



D. Transition densities 

Let us now look into individual modes of excitation. 
For this purpose we analyze the transition densities as- 
sociated with each excitation mode. We evaluate three 
kinds of transition densities [H, EH HE H3| 



(27) 



P?» = <*i|^(rWj(rl)|*o) = liV(f)i^(r), 



P t h q h (r) = ($,| i; q (r lty,(r f) l*o) = Yl M {r)P^ L {r) 



(28) 
(29) 
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FIG. 8: The same as FigJT] but for the B(IS2) strength func- 
tion in 20 O. 
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using the solution of the linear response equation ([21]) 
at an energy corresponding to a peak in the strength 
function. Here (? iq (r) is the usual particle- hole transition 
density while Pfq{r) and P^ h {r) are the particle-pair 
and hole-pair transition densities associated with the pair 
addition and removal amplitudes, respectively. 

Let us first discuss the giant dipole resonance (GDR) 
and the low-lying small peak seen in the B(El) strength 
function. We focus on 54 Ca, where we find two peaks 
at E x — 8.5 MeV and 16.0 MeV in the strength function 
shown in Fig[TJ (The corresponding peaks in the Landau- 
Migdal approximation are found at E x = 8.0 and 15.0 
MeV.) The transition densities evaluated at E x = 16.0 
MeV is shown in Fig(9]Ja). It is seen that the particle-hole 
transition density p p qL (r) is large around the nuclear sur- 
face, and that the neutron and proton amplitudes have 
the opposite phases. This is indeed the feature typical 
of the isovector giant dipole resonance (GDR). The char- 
acter of the excitation mode at E x — 8.5 MeV is dif- 
ferent from that of the giant resonance. This is seen in 
the transition densities shown in FigOJb), where we find 
a characteristic feature that the neutron amplitude of 
the particle-hole transition density dominates over the 
proton's in the exterior of the nucleus. The neutron 
amplitude has a node around the nuclear surface, and 
it exhibits significant magnitude also inside the surface, 
where the proton amplitude also has comparable magni- 
tude with the same phase as that of the neutron. This is 
the feature which is often interpreted as the soft dipole 
excitation or the pygmy dipole resonance characteristic 
to neutron-rich nuclei pj, El El III III ■ 

It is interesting to check how well the Landau-Migdal 
approximation can describe the transition densities and 
the mode characters. This is a non-trivial question since 
we have already seen that there is rather large differ- 
ence in the jB(E1) strength functions between in the full 
calculation and the Landau-Migdal approximation. The 
transition densities evaluated at E x =15.0 and 8.0 MeV 
are shown in the panels (c) and (d) in Figj9[ which are 
compared with those in (a) and (b). In both the cases of 
the giant dipole resonance ((a) vs. (c)) and of the soft 
dipole excitation ((b) vs. (d)), the basic features of the 
transition densities are the same in the full and LM cal- 
culations, although we see small but non-negligible differ- 
ences, e.g., in the relative size between the particle- hole 
transition density p ph (r) and the particle-pair transition 
density P pp (r) for the neutrons. This comparison sug- 
gests that the Landau-Migdal approximation can be used 
to describe the basic structure of these excitation modes 
while some reservation should be held when aiming at a 
quantitative description. 

Similar analysis of the transition densities is performed 
also for the isoscalar quadrupole modes in 54 Ca. We ana- 
lyze here the isoscalar giant quadrupole resonance having 
a broad peak around E x — 16.0 MeV and the low-lying 
quadrupole vibrational state peaked around E x = 2.5 
MeV. The transition densities evaluated at these peak 
energies are plotted in FigQI)] (a) and (b). The transi- 
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FIG. 9: The transition densities r 2 p^ L (r),r 2 P[^ L (r) and 
r 2 P^(r) for the isovector dipole response in 54 Ca. (a) and 
(b) are those in the full calculation, evaluated at E x = 16.0 
MeV (the peak energy of the GDR) and E x = 8.5 MeV (the 
peak energy of the soft dipole excitation), respectively, (c) 
and (d) are those in the Landau-Migdal approximation. We 
here plot the neutron amplitudes weighted with the volume 
element r 2 for the three kinds of transition densities. For the 
proton we plot only r 2 p? L (r). Note that we need to normal- 
ize the transition densities in terms of the _B(E1) strength of 
this modep^]. but we assume here the unit strength B(E1)=1 



tion densities of the corresponding peaks in the Landau- 
Migdal approximation are also plotted in (c) and (d). 
Comparing the full calculation and the Landau-Migdal 
approximation, we observe the same trends as found in 
the case of the giant dipole resonance and the soft dipole 
excitation. Namely there is no difference in the basic fea- 
tures of the modes, but quantitative details of the transi- 
tion densities depend on whether the velocity-dependent 
terms of the Skyrme interaction are taken account or not. 

We analyzed also the transition densities evaluated at 
the peaks of the dipole and quadrupole strength functions 
in 20 O. We obtained results similar to those in 54 Ca ex- 
cept for the low- lying peak at E x =8.5 MeV in the dipole 
response. The transition densities corresponding to this 
peak is different from those of the GDR nor those of the 
soft dipole excitation. We infer this peak as having non- 
collective nature. 
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FIG. 10: The same as Fig|9]but for the isoscalar quadrupole 
response in 54 Ca. (a) and (b) are those in the full calculation, 
evaluated at E x = 16.0 MeV (the peak energy of the ISGQR) 
and E x — 2.5 MeV (the peak energy of the low-lying surface 
vibration), respectively, (c) and (d) are those in the Landau- 
Migdal approximation. Note that we here normalized the 
transition densities in terms of the unit strength _B(IS2) = 1 
fm 4 . 



dependent terms by comparing with the reduced calcu- 
lation where the Landau-Migdal approximation is intro- 
duced to the residual interaction. The continuum QRPA 
using the Landau-Migdal approximation gives an overall 
correct description of the multipole responses, but influ- 
ences of the approximation are seen in the shift of the 
centroid energy by 2-3 MeV of the isovector giant reso- 
nances and also in the violation of the energy weighted 
sum rule by about 10-15 % in the case of SkM*. We 
found also small but non-negligible influence in the tran- 
sition densities for the low-lying dipole and quadrupole 
modes. Note however that the qualitative features of the 
transition densities are described well even in this case, 
and only very little influence is seen in the case of the gi- 
ant resonances. The analysis gives a partial justification 
for the use of the Landau-Migdal approximation. 
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IV. CONCLUSIONS 



DETAILS OF THE LINEAR RESPONSE 
EQUATION 



We have developed the continuum QRPA which is 
based on the Skyrme-Hartree-Fock-Bogoliubov energy 
functional. In deriving the residual interaction used in 
the QRPA, we have taken into account the velocity- 
dependent central terms (proportional to the t\ and t 2 
coefficients) of the Skyrme effective interaction in order 
to guarantee the energy weighted sum rule for the multi- 
pole responses, but we neglected the two-body spin-orbit, 
the spin-spin and the Coulomb interactions. 

The new continuum QRPA is applied to the isovector 
dipole and the isoscalar /isovector quadrupole responses 
of medium-mass neutron rich nuclides 20 O and 54 Ca us- 
ing the SkM* parameter set. It is confirmed numerically 
that the energy weighted sum rule is satisfied up to the 
enhancement factor relevant for the isovector responses. 
This is because the velocity dependent terms are taken 
into account in a way consistent with the Hartree-Fock- 
Bogoliubov mean-fields of the ground state. We thus 
constructed the first Skyrme continuum QRPA formal- 
ism that satisfies the sum rule. 

We have also examined the importance of the velocity 



Here we give some details of the linear response equa- 
tion (|2"T]) . The functions a qq i , b, 



qq' ! vqq' ; Cqq' 



and d q appearing 

in Table U and Eqs.© and (TIT))) are expressed as follows 
in terms of the effective interaction parameters and the 
local densities: 



\t (l - Xq) 

+(a + 2){a + l)^*s(l + \^)p a 

-T2h(x3 + \) 

x a(a-l)p Q - 2 E 9 -P'" 

+4ap a - 1 p q + 2p a 

*0(1 + \x ) 

+(a + 2)(q + l)^t 3 (l + \^)p a 
x (a(a - \)p°- 2 V „ pi, + 2ap a 



lqq> 



= < 



(? = </') 



b qq > = 



(q ^ q') 

-^(ti(i-ai) + 3ta(H-x2)) (q = q') 

-I{i x (l + \xi)+t 2 {l + \x 2 )} {q^q') 
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Spa 



2 E_0i_Sf [^S PlL /r 2 + S po f qL 



[Sp q ] L 

[Sp+q]L 
[Sp-q]L 



2x T^"ftE ' b qq '[ 5 P q '\L 



n 

2 m 



2 x -ftT-AzE,' V'[<W]i 



2x ir |p?E,' a qq ,[8p ql ] L + p q a q [8p + , q ] L + 2{A Pq - 2r 9 ) ^ fwP/vli 

c „„o (.f^Zisvp^}^- 1 - ,f^[SVp q 



'P?']i 



2L+ 



[•JVpj^^ 1 2 x 2=£ { V fovk + p* ( V - g£r<w M^]^" 1 + ^Pp, <w L+1 } 

[SV Pq ]l= L+1 2 x ^ {-y^^ J2 g , b qq ,{5p q ,] L + ^^. Pq j2 q , c qq ,[8Vp ql ]>r L -'+p q - ^c^o^v^]^ 1 } 

i^iijr^ 1 -2 x ^ {s^xp.e^ vi^jvi^- 1 - 4SPp«e,/ v[<52y>r i+1 } 

[St,]l 2 x ^ {r,^, b qq ,[8 Pq ,] L + ±^£ ? ,(V - <w) (y/^fr[SVp*]i= L - 1 ~ P*]^ 1 ) } 

TABLE II: The second term of r.h.s. of the linear response equation (|21[) 

f -±(t 1 {x 1 -l)+9t 2 (x 2 + l)) (q = q') - Vb 

I + - 3t 2 (l + 5Z2)} (9 ^ -?') 
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